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Abstract. Computer simulations of inhomogeneous soft matter systems often require 
accurate methods for computing the local pressure. We present a simple derivation, 
based on the virial relation, of two equivalent expressions for the local (atomistic) 
pressure in a molecular dynamics simulation. One of these expressions, previously 
derived by other authors via a different route, involves summation over interactions 
between particles within the region of interest; the other involves summation over 
interactions across the boundary of the region of interest. We illustrate our derivation 
using simulations of a simple osmotic system; both expressions produce accurate results 
even when the region of interest over which the pressure is measured is very small. 



PACS numbers: 62.50.-p,83.10.Rs,82.39.Wj 
1. Introduction 

Molecular dynamics (MD) simulations are widely used to study spatially inhomogeneous 
soft matter systems. In such simulations, it is often necessary to compute the local 
pressure in a small region of the simulation box, containing only a few atoms or 
molecules. Examples include calculations of interfacial free energies [I], measurements 
of osmotic pressure gradients [2j [3] and tests of coarse-grained hydrodynamic theories 
[3]. While accurate expressions for the local pressure exist, their derivation is rather 
involved. In this paper, we present a much simpler derivation, which leads to two 
equivalent expressions for the local pressure. One of these expressions is analogous to 
the local stress tensor of Lutsko [5] and Cormier et al. [6]; the other is, to our knowledge, 
new, but is similar in spirit to the "Method of Planes" of Irving and Kirkwood [7] and 
Todd et al. [Sj. We show that both these expressions give accurate results for the local 
pressure in soft matter systems, even when computed over very small spatial regions. 
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2. Background 

For a spatially homogeneous, closed system, the pressure is commonly computed by 
taking the average of an instantaneous "pressure function" : 

where N, V and T are the number of particles, volume and temperature, ks is 
Boltzmann's constant, r*j and fj are the positions of particles i and j, fy = fi — r}, 
and fij denotes the force exerted on particle i by particle j [9j [TO]; the double sum 
runs over all pairs of particles, avoiding double counting. The first and second terms in 
Eq. (|T]) arise from the kinetic energy of the particles and from interparticle interactions, 
respectively. Expression ([I]) can be derived in a few steps starting from the Clausius 
virial relation 

(e^+X>^) = o, (2) 

\i=l "H j=l / 

in which pi and are the momentum and mass of particle i, and Fi is the total force 
acting on particle i, due to other particles and the walls of the container [HI IT2"] . 

For spatially inhomogeneous systems, one can measure directly the pressure across a 
local plane within the simulation box via the "method of planes" (MOP), first proposed 
by Irving and Kirkwood [7] and later rederived by Todd et al. [8] ; this works well when 
the plane over which the pressure is computed is large, but leads to poor statistical 
sampling when computing the local pressure in a small region (i.e. over a small plane). 
Alternatively, the pressure in a local region of interest can be measured using a local 
version of Eq.((T]). This has the advantage that the region of interest can (in principle) 
be of arbitrary shape, and that statistical averages are taken over a volume rather than 
an area. The following local pressure function was proposed by Lutsko [5] and later 
reformulated by Cormier et al. [6] (note that these authors considered the full stress 
tensor, while, for simplicity, we focus only on the scalar pressure, which we assume to 
be locally isotropic): 

1 / N \v\ 2 N ~ x - \ 
P(f) = ^ ( E + £ ■ r^h ■ (3) 

\i=l 1 i=l j>i I 

Here, Q is the volume of the region of interest, centred on r, Aj is unity if particle i lies 
within the volume fi, and zero otherwise, and Uj is the fraction (0 < 1^ < 1) of the line 
joining particles i and j that lies within Q (see Figured]). The local pressure expression 
(J3]) is analogous to the global one ([1]), but with two important differences. Firstly, in 
the kinetic part of Eq.OH]), only those particles which are inside the region of interest (at 
time t) are included. Secondly, the components of the interaction term are weighted by 
the fraction of the line joining particles i and j that is inside the region of interest; this 
highlights the crucial importance of correctly accounting for interparticle interactions 
which cross the boundary of the region of interest. Note that particles i and j may both 
be outside the region yet still contribute to the interaction pressure; see Figure CD 
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Figure 1: Contributions to the interaction part of Eq.([3]) from different pairs of particles: 
the fraction Uj of the line joining particles i and j which lies inside the region of interest 
is shown as a solid line, bij denotes the position where the line joining particles % and j 
crosses the boundary. 

Both Lutsko [5] and Cormier et al. [S] derive the stress tensor version of Eq.([2]) 
by Fourier transforming the continuity equation for the local momentum flux and 
applying Newton's second law. In this paper, we present a simpler and arguably more 
intuitive, real space derivation, which is directly analogous to the derivation of the 
global pressure expression, Eq.flT]), from the Clausius virial relation (j2J). This derivation 
leads both to Eq. ([3]) and to a new expression for the local pressure involving the flux of 
particle momentum across the boundaries of the region of interest, together with cross- 
boundary interactions. The equivalence between these two expressions shows directly 
how the relation between surface flux and volume pressure measurements extends to the 
atomic scale. Our approach may also prove useful in future for deriving local pressure 
expressions for systems with different dynamical rules, such as run-and-tumble swimmers 
[T3] or particles with viscous dynamics [10] . 

3. Derivation of expressions for the local pressure 

To compute the local pressure at some position r in the system, we consider a local 
region of fluid, of volume Q, centred on r. Particles within this region interact with other 
particles both inside and outside the region. During a given time interval, particles will 
enter and leave the region of interest. 

3.1. The Schweitz virial relation 

We begin with an analogue of the Clausius virial relation fl2]), derived by Schweitz, for 
open systems [H]. The Schweitz virial relation states that 
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where, as above, the function Aj(t) measures whether or not particle i is within the 
region of interest at time t, and its time derivative dAi(t)/dt = Aj produces a positive 
or negative 5-function peak at the moment when particle % enters or leaves the region 
of interest 0. The first term in Eq.(j3J), fetLi ^r-A-j(r)) = ^Un, is the average kinetic 
energy of particles in the region of interest and is directly analogous to the first term 
in Eq.([2J). The second, interaction, term is analogous to the second term in Eq.(j2]) - 
we assume that there are no external forces so the total force on particle i is given by 
the sum of interactions with all other particles in the system. The final term in Eq. (pE}, 
J^iLii^i -Pi)Aj\ = $, accounts for the exchange of particles between the region of 
interest and its surroundings. Particles entering the region contribute r*j • while those 
leaving contribute —r\ ■ pf, these do not cancel because the momentum vectors Pi for 
particles entering and leaving are opposite in sign. 

We now split the second term in Eq.fjlJ) into contributions due to interactions with 
particles inside and outside the region of interest: 

X» (EAj) A i (f) = V iBt + V« Bt , (5) 
i=i \j# J 

where V int = (Ef=iE^ 4 n ■ fa AiAj) = (E 1 =4 1 Ej^ij ■ fa A*Aj) = (Einin^i ' fap 
contains contributions where both particles are inside the region of interest and V ex t = 
(EfcLiEjyi^ " fa Aj(l — A.j)\ = ^Einout^i ' fij) contains contributions where particle i 
is inside the region and particle j is outside. Substituting Eq.flS]) into Eq.pj allows us 
to write the Schweitz virial relation as 

E kin + V int + Vext + $ = 0. (6) 



3.2. Expressions for the local pressure 

We next use the Schweitz virial relation to derive expressions for the local pressure. The 
local pressure has two components: a kinetic component, Pkinif), which is given by the 
normal flux of particle momentum across the boundaries of the region of interest and an 
interaction component, Pintif), which is the surface density of interparticle forces across 
the boundary. Throughout this paper the normal to the boundary is assumed to point 
in the outward direction. 

X For closed systems, the Clausius virial relation, Eq.([2|), can be derived by setting the time derivative 
of the function G = ^E^Li^'P*^ to zero in steady state. For open systems, the derivation 
of the Schweitz virial relation follows a similar route, but takes into account the contributions to 
Giocai(r) = (j2iLi(^i ' Pi)^i{'fj) from particles entering and leaving the system [14] . 
§ here we have used the fact that EiliEj^i^ ' fa = YaJi Y.jyi^j - fij, since fij = ~fji 
[ID] . We have also introduced a new notation: EininS'y = Eili 1 Ej>i A i A jSij- Likewise, we 
define Eoutout&j = EiL^EixC 1 ~ A «)(l - ^j)dij and EinoutSy = E^Ii E^i A *(! - A j).9y = 
El'i 1 EixIMl - A^gtj + A,(l - A*)^]. 
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Kinetic component The kinetic component of the local pressure can be related to the 
component $ = (j^iLi " Pi) ^i) of the Schweitz virial relation. To see this, we split 
the particle momentum, pi, into its components normal and tangential to the boundary: 
Pi = {Pi ■ n)h + (pi ■ to give 



N 



$ = ( J2 A i (Pi ■ n)(fi ■ n) + ■ t)(fi -t) ). (7) 



=1 



Assuming that the density of particles is uniform across the region of interest, the second 
term in Eq.fjTj) averages to zero. Next, we notice that because A, is nonzero only when 
a particle is at the boundary, r*j • n may be taken outside the summation. Assuming 
(without loss of generality) that the region of interest is cubic with the origin at its 
centre and sides of length L, fi-n = L/2, and the total (outward) momentum flux 
(j2i(Pi " ^)Aj^> across each of the 6 faces is —Pki n L 2 . Eq.flT]) therefore reduces to 

$ = -3nP kin (r). (8) 

Interaction component In a similar way, the interaction component of the local 
pressure, Pmt(r), can be related to the component V ex t of the Schweitz virial relation. 
We first note that the position vector r*j may be written as f\ = fey + fijhj, where fe^ 
denotes the position where the line linking particles i and j crosses the boundary of the 
region of interest, and 1^ is the fraction of the line linking particles i and j that lies 
inside the region of interest (see Figure [1]). V ex t is then given by 

Vext = ( £ h ■ fij + kj (nj ■ fij) ) (9) 
\ in out / 

E (fij ■ h ) ■ fl )) + c t + \T, l a ( r a ■ fa. 



\in out 



where the second line follows from splitting the interparticle force / y - into its components 
normal and tangential to the boundary: = (fy ■ h)h + (fy ■ t)t, and C t = 

(j2inout(fij " t) i^ij 't)\ Focusing on the first term, we note that points to the 
boundary and so (assuming the same cubic geometry as above), fo y - • h = L/2. Denoting 
as a the average outward normal force per unit area crossing the boundary, due to the 
(zin, jout) interactions, we obtain 

/ E (fiJ ■ n) (kj ■ n) ) = -6 (%) (L 2 a) = -3Qa. (10) 



\inout / 

The interaction component of the pressure, P int (r), is equal to the total surface density 
of normal force crossing the boundary: Pi nt = a + £, where £ is the normal contribution 
due to pairs of particles i and j which are both outside the region of interest (see Figure 
[Q. We demonstrate in the Appendix that 



3Cl£ = C t -( E l ij(rij-fij)) (li; 
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so that, putting together Eqs(19])- (fTT]) . we can relate P int to V ext by: 

- 30P int (fO = - 3Q(a + = V ext - V corr} (12) 

where 

V corr = ( E ( r *i ' fij) E ifv ' A?) / ' (^) 



\ in out out out 



"Boundary" expression for the local pressure An expression for the local pressure, P(f) , 
can be obtained by combining Eqs.flH} and f|T2|) : 



P(r) = P hi n + Pint = ~ [$ + V ex , - V c 



1 

3fi 




'"i ' Pi) Aj -|- ^ ' (rj ^ij^ij) ' fij ^ ] liil 
in out out out 




Eq.( IT4|) provides a simple prescription for computing the local pressure. The first term 
sums over all particles which enter or leave the region of interest and is equivalent 
to the momentum flux density due to particles crossing the boundary, while the 
remaining terms, which account for the force density at the boundary due to interparticle 
interactions, sum over all pairs of particles for which the line connecting the two particles 
crosses the boundary of the region of interest. In the case where the region of interest is 
large, P(r) is dominated by the contributions of <E> and V ext {V corr becomes negligible); 
however, as we show below in Figure El V CO rr makes an important contribution when the 
region of interest is small. Eq.( jT4l) provides an alternative to existing local pressure 
expressions, and demonstrates explicitly how the relation between surface flux and 
volume pressure expressions extends to very small regions of space. 



"Volume" expression for the local pressure 

The Schweitz virial relation provides a direct route from this "boundary" expression to 
the more usual expression for the local pressure, which involves a sum over particles 
within the region of interest. Inserting Eq. (TH|) into the Schweitz relation (JBJ), we obtain: 

P{r) = ^ [E kin + V int + V corr ] (15) 

= ok (e ^r A *( f ) + E ? a -fn + E h (4- ■ h) + E h ■ fa) ) 

°" \i=l "h inin in out out out / 

-4 (e^a.M + ££^ ■/«)), 

\i=i 1 i j>i i 

where the last line follows from the fact that Uj = 1 if both particles are inside the 
region of interest (assuming the boundary is everywhere concave). Eq.f ll5p is identical 
to Eq.Q, and constitutes a local version of the global instantaneous pressure function, 
Eq.flTJ). 
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4. Molecular dynamics simulations 

We now illustrate, using molecular dynamics simulations, the Schweitz virial relation 
(jUJ) as well as the two expressions for the local pressure, Eqs. ffT4"l) and ffT5j) . In our 
simulations, a periodic box contains 5000 particles at a density p = 0.8 (reduced 
units [IQ]), which interact via a repulsive Weeks, Chandler, Andersen (WCA) potential 
[15] (a truncated and shifted Lennard- Jones potential). The particle size a = 1, the 
interaction parameter e = 1 (both in reduced units) and the box length is 18.42cr. The 
system is simulated using the velocity Verlet algorithm [10] with timestep At = 0.001 
(reduced units) and is maintained in the canonical (NVT) ensemble using a Nose- Hoover 
thermostat at temperature, T = 1.0 (reduced units). All runs are equilibrated for 4 x 10 4 
timesteps prior to data collection, and data is collected over at least 2 x 10 6 timesteps. 
Errors are computed using bootstrapping [16J, with 1000 replica datasets. 



4-1. A homogeneous fluid 

We first consider a homogeneous fluid, for which the local pressure, measured using 
Eqs.( !T4"]) and (JT5]) . should match the global pressure, measured using Eq.([T]). We define 
a cubic region of interest, located in the centre of our simulation box, whose size L we 
vary from L = 5.5a to L = 3a. For this smallest value of L, the region of interest 
contains only ~ 21 particles on average. 
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Figure 2: Molecular dynamics simulation results for the homogeneous fluid, (a): 
Ekin + Vint + Vext + $ as a function of the size L of the region of interest, (b): Local 
pressure in the region of interest, computed using Eq. (TT4"|) (blue squares) and Eq. (fT5l) 
(black circles). The apparent correlation between these results arises because the same 
simulation data set was used in both cases. The green line shows the global pressure 
computed using Eq.([T]). The inset shows the pressure computed using Eq. (}T5]) . not 
including V corr (red triangles; note the very different scale on the pressure axis). 



Figure |2k shows E kin + Vi nt + V ext + $ as a function of the size L of the region of 
interest. As predicted by the Schweitz virial relation ()6]), this quantity is zero within 
our error bars (note that of the individual terms in this sum, two are ~ 5 and the 
other two are ~ 0.8). Figure [2b shows the local pressure P(r) in the region of interest, 
computed using expressions ( Til"]) and ( IT5l) . as a function of L. Both expressions give 
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results in excellent agreement with the global pressure across the whole simulation box, 
using Eq. ([1]). The inset to Figure [2b, which shows Eq.f lT5|) . neglecting the V corr term, 
demonstrates the importance of correctly accounting for the boundary terms: neglecting 
V corr gives a large error, which increases as L decreases. 



4-2. An osmotic system 
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Figure 3: Molecular Dynamics simulations for the osmotic system, (a): E^ n + V% n t + 
Vext + as a function of the concentration N u /V u of solute particles, computed for a 
small subregion (of dimensions L x = L y = L z = 6a) within the solution region (black 
circles) and for a small subregion (L x = 2a, L y = L z = 7a) outside the solution (blue 
squares), (b): Osmotic pressure difference II between the solution and its surroundings, 
computed using Eq. (fT4l) (green triangles) and using Eq. (fl5l) (blue circles) (for the same 
subregions). Results are also shown for a direct computation of the normal solute- 
membrane force per unit area (red squares). Error bars are smaller than the symbols. 
The prediction of the van't Hoff equation, II = (N u IcbT) /V u , is also shown (black line). 



We next consider an osmotic system, in which the local pressure differs in different 
parts of the simulation box. To construct this system we label a subset N u of the 
particles "solute" and the remaining particles "solvent". Solute particles are confined 
to a cubic region of volume V u ~ 786er 3 in the center of the simulation box by a smooth 
confining potential; solvent particles do not experience this potential and are free to 
move throughout the simulation box. The confining potential acts as a semi-permeable 
membrane, resulting in an osmotic pressure difference, IT, between the "solution" region 
where the solutes are confined and the rest of the simulation box. 

To compute the osmotic pressure, we define local regions of interest inside and 
outside the solution compartment (for dimensions see the caption of Figure [3]). Figure 
[3^ shows that the Schweitz virial relation fl6]) is obeyed in both of these local regions, 
over a range of solute concentrations. The pressure in the two local regions can be 
computed using Eqs. f[T4|) and f[T5l) ; subtracting the results gives the osmotic pressure 
difference, II. Figure [3b shows II, as a function of the concentration N u /V u of solute 
particles in the solution compartment; both methods for computing the local pressure 
produce results in excellent agreement with a direct calculation of the osmotic pressure 



Computing local pressure 



9 



obtained by measuring the average confining force on the solutes, per unit area of the 
confining box. 

5. Conclusions 

Accurate methods for computing the local pressure are essential for simulating 
inhomogeneous soft matter systems. In this paper, we have derived two equivalent 
expressions for the local pressure in a molecular dynamics simulation. The "boundary" 
expression, Eq.( TH|) . is, to our knowledge, new. This expression involves summation 
over interactions between particles within and outside the local region of interest and 
is similar in spirit to the "method of planes" approach of Irving and Kirkwood [7] 
and of Todd et al. [8J. The "volume" expression Eq. f|T5|) is a local analogue of the 
function commonly used to compute the global pressure in homogeneous simulations; 
this involves summation over interactions between pairs of particles within the region 
of interest. This expression was previously derived via a Fourier transform method by 
Lutsko [5 J and by Cormier et al. [6] ; our derivation, based on the Schweitz virial relation, 
provides a simple real-space alternative. Importantly, both local pressure expressions 
take careful account of interactions close to the boundary: this is crucial when the region 
of interest is of the order of the particle size. 

The derivation presented in this paper, via the Schweitz virial relation, 
demonstrates explicitly how the equivalence between surface flux and volume expressions 
for the pressure, familiar from macroscopic systems, plays out on very small (atomistic) 
lengthscales. This approach may also prove useful in future for deriving local pressure 
expressions in systems whose dynamics are more complex: for example systems with 
viscous dynamics [10], or active matter systems in which particles are self-propelled 
and/or chemotactic [13] . Here, the Fourier transform method of Lutsko and Cormier et 
al. [3, [6] might prove challenging, but we hope that our real-space method should hold 
with only minor modifications. 

Finally, we note that an important assumption made in this work is that the local 
pressure is isotropic: we therefore derive expressions for the scalar pressure rather than 
the local pressure tensor, as in previous work [5], EE IB] ■ We believe that it should be 
possible to extend the present derivation to obtain the analogous expressions for the 
pressure tensor, via a tensor analogue of the Schweitz virial relation. For the present, 
we leave this as an interesting avenue for future work. 
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Appendix: Pairs of particles which are both outside the region of interest 



Here, we discuss the contribution to the interaction component of the pressure made 
by pairs of particles i and j, in the cases where both particles are outside the region of 
interest, but (as illustrated in Figured]), the line joining particles i and j crosses the 
boundary. Our aim is to derive Eq. ( TTTj) . Assuming that the region of interest is cubic 
(with side length L and origin at the centre), the line joining particles i and j crosses 
the boundary on two different faces. We define these crossing points by the position 
M nr ,A T?;n-nT-o PTk iiinD+rotoo fV, of 7 — 7-i( 2 ) — TM> This allows us to write 



vectors b\j' and bU . Figure [Jji illustrates that l^r. 



V - b ( 

*3 V 



\out out 
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We now resolve b\j and b\j into components normal and tangential to the boundary 



on faces 1 and 2 respectively [b. 



0® 



n 



(1 ))nW + (W ■ i^)tW etc]. Noting that 



« . aw 



(6 



*(2) 



n 



(2)> 



L/2, we obtain 




3f2f - (A (1) + A 1 



(•2) 



\out out 

where £ is the (outward) normal force per unit area crossing the boundary due to particle 
pairs which are both outside the region of interest, 

i(2) 



and D\ z> = (Eoutout(&lf - i {2) ){fji - i {2) )) ( we use = "A?)- Assuming that particles 
are homogeneously distributed across the region of interest, D t (1) = = D t /2 (by 
symmetry, see Figured^). 
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Figure 1: (a) Geometric illustration of the case where particles i and j are both outside 
the region of interest, (b) Illustration of the simplified one dimensional scenario used to 
show that C t = —D t . 



We now demonstrate that D t = —C t , where C t , as defined in Section [3] is the 
tangential component of the (i in, j out) contribution to V ext - We study the simplified one 
dimensional scenario shown in Figure[Tb, in which we consider only particles i and j lying 
along two lines parallel to the x axis. We further assume that interactions are restricted 
in range, such that each particle % interacts with only two interaction partners j, as 
shown in Figure [lb, with tangential forces / in opposite directions. These interactions 
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contribute respectively f(xi — (3) and —f(xi + (3) to the sum YU ' t)(fij ' We 

now sum the contribution of the interactions shown in Figure [lb for the cases where 
particle i is inside the region of interest (— L/2 < X{ < L/2) and the line joining i and 
j crosses the top face of the region - i.e. the contributions to C t for this face. The 
result is f/3d(/3 — L), where d is the number of particles i per unit length along the line 
[Jjl Next, we sum the contributions for the cases where particle i is outside the region 
of interest, but the line joining i and j still crosses the top face of the region (i.e. the 
contributions to D t for this face). Contributions ±(x =F 0)f are made by particles i 
for which {-L/2 - (3) < Xi < -L/2 or L/2 < Xi < (L/2 + 0) (case C in Figure Eb). 
Integrating over these ranges of Xi, we obtain a total contribution to D t of f/3d(L — (3). 
Thus the contributions to the (i in, j out) tangential term C t are exactly compensated by 
the contributions to the (zout, j out) tangential term D t , and Eq.OI) reduces to Eq. ( TTTTl . 
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|| This result follows from noting that particles i for which {—L/2 + ft < xi < L/2 — ft) have 2 partners 
and contribute —2fft (case A in Figure[T}D), particles for which —L/2 < X{ < — L/2 + ft have 1 partner 
and contribute — (xi + ft)f while particles for which L/2 — ft < xt < L/2 have 1 partner and contribute 
(xi — ft)f (case B in Figure [TJd); these contributions are then integrated over —L/2 < x < L/2: 

df 



- j_^ +P (x + ft)dx- n; 2 ^ 2ftdx + tij 2 % (x -ft)dx = fftd{ft - L) . 



